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The planar bistable device [Tsakonas et al, Appl. Phys. Lett., 2007, 90, 111913] is known to have two distinct 
classes of stable equilibria: the diagonal and rotated solutions. We model this device within the two-dimensional 
Landau-de Gennes theory, with a surface potential and without any external fields. We systematically compute a 
special class of transition pathways, referred to as minimum energy pathways, between the stable equilibria that 
provide new information about how the equilibria are connected in the Landau-de Gennes free energy landscape. 
These transition pathways exhibit an intermediate transition state, which is a saddle point of the Landau-de 
Gennes free energy. We numerically compute the structural details of the transition states, the optimal transition 
pathways and the free energy barriers between the equilibria, as a function of the surface anchoring strength. 
For strong anchoring, the transition pathways are mediated by defects whereas we get defect-free transition 
pathways for moderate and weak anchoring. In the weak anchoring limit, we recover a cusp catastrophe situation 
for which the rotated state acts as a transition state connecting two different diagonal states. 


I. INTRODUCTION 

Nematic liquid crystals are complex anisotropic liquids 
with long-range orientational ordering |[l|]. Nematics in 
confinement present a whole host of new theoretical and 
applications-oriented questions focussed on the complex 
inter-relationship between material properties, geometry, tem¬ 
perature, boundary effects and external fields in equilibrium 
and non-equilibrium phenomena. From a modelling point of 
view, we have new and challenging questions on pattern for¬ 
mation, defects, interfacial phenomena and rheology. From a 
technological point of view, nematics in micron-scale or nano¬ 
scale geometries open new doors for optical and display ap¬ 
plications, and more recently novel biological insight y,[j]. 

We re-visit the planar bistable device first reported by 
Tsakonas and his co-workers il. The planar bistable device 
has been well-studied in the liquid crystal community in re¬ 
cent years iil] , partly because of its geometrical simplicity 
and partly because of its rich modelling landscape. It typi¬ 
cally comprises a periodic array of shallow square or rectan¬ 
gular wells filled with nematic liquid crystalline material. The 
well surfaces are treated to induce tangent or planar boundary 
conditions so that the nematic molecules, in contact with the 
well surfaces, are constrained to be tangent to the well sur¬ 
faces. Tsakonas et al. [0] adopt a two-dimensional Landau-de 
Gennes modelling approach and focus on planar nematic equi¬ 
libria on the bottom square or rectangular well cross-section, 
on the grounds that non-planar three-dimensional equilib¬ 
ria are energetically expensive for shallow three-dimensional 
wells and hence, not physically relevant. They compute local 
minimizers of a two-dimensional Landau-de Gennes energy 
on a square or a rectangle with tangent boundary conditions 
and obtain two distinct classes of optically contrasting ne¬ 
matic equilibria: the diagonal solutions and the rotated solu- 
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tions and their modelling results are supported by parallel ex¬ 
perimental work. The tangent boundary conditions naturally 
create a mismatch in molecular alignments or defects at the 
square vertices. As the name suggests, the nematic molecules 
roughly align along the square diagonal for the diagonal so¬ 
lution. The rotated solutions have a more distorted profile as 
simulations and optical data suggest that the molecules rotate 
by about 180 degrees between a pair of parallel square edges. 
There are two diagonal solutions, one for each square diago¬ 
nal, and four rotated solutions, related to each other by a 90 
degree rotation. Luo et al. (H] build on the work in |0 and 
carefully study the dependence of the diagonal and rotated so¬ 
lutions on the surface anchoring strength, denoted by a surface 
anchoring coefficient W. In particular, they find that the ro¬ 
tated solutions only exist above a certain critical anchoring 
strength and the critical anchoring strongly depends on the 
material parameters and temperature. 

Here we take the modelling work further by a systematic 
numerical investigation of the free energy pathways of the 
planar bistable device within the framework of the Landau- 
deGennes theory. This is a significant forward step since we 
not only recover the free energy minimizing states but also 
compute structural and energetic information about the tran¬ 
sient states connecting pairs of distinct free energy minima. 
As we vary the surface anchoring strength, we find three dis¬ 
tinct regimes: the strong anchoring, moderate anchoring and 
weak anchoring regimes. In the strong anchoring limit, we re¬ 
cover the familiar diagonal and rotated solutions as local free 
energy minimizers. As noted in 10,0, there are two differ¬ 
ent diagonal solutions and four distinct rotated solutions. We 
compute the transition pathways between the different stable 
solutions. In particular, we focus on the so-called minimum 
energy pathways, where every point on such a pathway is an 
energy minimum in all but a single distinguished direction in 
the phase space. There are a number of interesting findings to 
highlight here. Firstly, each minimum energy pathway is fea¬ 
tured by a transition state, which is a saddle-point, an unstable 
critical point of the Landau-de Gennes energy. We identify 
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the transition states as being local free energy maxima along a 
minimum energy pathway connecting free energy minima. In 
the strong anchoring limit, the transition states exhibit ±1/2 
defects along the square edges. There may be multiple mini¬ 
mum energy pathways between a pair of distinct minima and 
these transition pathways can contain different number of de¬ 
fects. We interpret the optimal transition pathway as being the 
minimum energy pathway with the smallest free energy bar¬ 
rier and the least energetically expensive transition state. It is 
noteworthy that we do not find a direct minimum energy path¬ 
way between pairs of diagonal solutions or pairs of rotated 
solutions and all minimum energy pathways connect diagonal 
and rotated solutions. In the moderate anchoring regime, the 
transition pathways are not mediated by defects but are rather 
mediated by localized anchoring breaking along the edges, 
which induces a global transition between diagonal and ro¬ 
tated solutions. In the weak anchoring regime, the rotated so¬ 
lutions cease to be locally stable (consistent with the numeri¬ 
cal findings in Jst]) but do exist as transition states connecting 
the di^onal solutions. In particular, the numerical methods 
in 10, 0 cannot capture the persistence of rotated solutions, 
as non energy-minimizing solutions, in the weak anchoring 
regime. 

We make some remarks on the novelty and importance of 
our approach. The transition states are not stable but they may 
be experimentally observable. Indeed the transition states in 
the strong anchoring limit are reminiscent of recent experi¬ 
mental observations of metastable states with point defects 
in shallow nematic chambers Jsl] . Further, our numerical re¬ 
sults demonstrate a subtle dependence of the optimal transi¬ 
tion pathways on the anchoring strength. In the strong anchor¬ 
ing limit, the optimal transition pathway is almost indepen¬ 
dent of the anchoring strength whereas the free energy barrier 
decreases monotonically with anchoring strength in the mod¬ 
erate and weak anchoring regimes. Switching mechanisms 
rely on a complex interplay between anchoring strength, ma¬ 
terial properties and external fields and one could use our nu¬ 
merical methods to compute and analyze optimal switching 
pathways in the presence of external fields. 

The paper is organized as follows. In Section |II1 we re¬ 
view the Landau-de Gennes theory for nematic liquid crys¬ 
tals. In Section nni we compute the Landau-de Gennes free 
energy minimizers and in Section [IVl we compute the tran¬ 
sition states and the transition pathways between free energy 
minima, focussing on three different anchoring regimes. We 
conclude in Section IVl with future perspectives. The compu¬ 
tational methodologies used in this paper are described in the 
appendix. 


be a square in the (x, -plane, which defines the bottom sur¬ 
face of the well, 

{ix,y) :0 <x,y <L} , (1) 

where L denotes the system size. As in Jsl], for two dimen- 
sions, the Q tensor can be represented by a symmetric trace¬ 
less 2x2 matrix. 


Q = 


Qll Qi2 
Qi 2 —Qll 


s{2n (g) n — I), 


( 2 ) 


where n is the director or the distinguished direction of molec¬ 
ular alignment and 5 is a scalar order parameter that measures 
the degree of orientational ordering about n. 

We work with a simple form of the Landau-de Gennes en¬ 
ergy as given below, with no external fields. 


Tr = 

± 


JQ 


-aTrQ2 - yTrQ^ 2(TrQ2)2 


dA (3) 


/ ^|VQ|2dA+ [ W^|(Qii,Qi2)-qpda. 
Jn ^ JdA 


The first integral is the bulk free energy that drives the 
nematic-isotropic transition as a function of the temperature 
Cllil. The coefficient a is the re-scaled temperature; we 
work with temperatures below the critical nematic supercool¬ 
ing temperature and hence, a > 0; the coefficients B^C are 
positive material-dependent constants. For a two-dimensional 
Q as in (|2]), trQ^ = 0, and the bulk free energy simplifies 
to the familiar Ginzburg-Landau potential iH]. The second 
term is the one-constant elastic energy and Hei > 0 is an elas¬ 
tic constant; there are more general quadratic elastic energy 
densities but we believe that the one-constant energy density 
suffices for qualitative purposes. The third term is a Durand- 
Nobili surface anchoring energy El that enforces the pre¬ 
ferred tangential anchoring on the square edges, with anchor¬ 
ing strength W, reference configuration q and boundary el¬ 
ement da. There are multiple choices for the surface energy 
(see (I] for comparisons between three potential candidates 
for the tangent surface energy) but the Durand-Nobili energy 
has the desired numerical stability properties for all relevant 
ranges of W. 

To reduce the number of input parameters in the model, we 
now rewrite the free energy functional in its dimensionless 
form. Substituting Eq. ^ into ^ and defining x = xjL, 
§ = kei = nei/aL‘^, = CQVa, W = 

WjaL, we obtain 


II. COMPUTATIONAL MODEL 

Following the paradigm in li a, we model the pla¬ 
nar bistable device within the two-dimensional Landau-de 
Gennes theory. In Landau-de Gennes theory, all the micro¬ 
scopic details about molecular shape and interactions are aver¬ 
aged out, and the nematic state is described by an macroscopic 
order parameter, Q. We take the computational domain, fl, to 


^ = / (Qn + Qi 2 - IrdA 

/a 


= [ (Q! 

Ja 

L 

[ WKQii, 

JdA 


|VQii|2 + |VQi2| 


dA 


(4) 


+ 


Q 12 ) - q\^dd. 


Here, there are only two input parameters: the dimensionless 
elastic constant kei and the dimensionless surface anchoring 
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strength W. We drop the tildes in the subsequent text and 
all results are to be interpreted in terms of the dimensionless 
variables. 

We need to prescribe a suitable form for the reference con¬ 
figuration, q, that enforces the tangent boundary conditions 
on the edges. Following the formulation in Luo et al. (I], 
wesetq(x,^) = 5(x, (cos 26>(x, ^), sin 26>(x, ^)), where 
= '7r/2, 0{x,0) = 0{x,l) = 0, and 
s(t,0) = s(t, 1) = s(0,t) = s(l,t) = /(t). The scalar s 
has to vanish at the vertices with strong tangent anchoring, 
because of the mismatch in 0 at the vertices. Hence, we define 


fit) = < 


t/d, 

1 , 

(1 — t)/d 


0 < t < d, 
d < t < 1 — d, 
l-d<t<l, 


with d chosen to be In the liquid crystal literature, 

is known to be typically proportional to defect core sizes 
illfl^ and hence, we take / to be approximately unity every¬ 
where (a minimum of the bulk free energy) except for small 
neighbourhoods around the square vertices, where we expect 
to see defects. 

We can use standard methods from calculus of variations to 
derive the Neumann boundary conditions on the edges (from 
the Durand-Nobili surface energy in Eq. dH)), 


^Qcx(3 

dx 

dx 

dQaj3 

dy 

dQap 

dy 


W 

- {Qa /3 - qa/d) on X = 1, 

^el 

W 

H- iQap - Qap) on a; = 0, 

^el 

w 

- (Qa/3 - qap) on y = 1, 

l^el 

w 

H-(Qa/3 - qa/3) on y = 0. 

I^el 


We use Greek indices to denote the components of the Q- 
tensor, = Qu and Q 12 above. 

The numerical methods for the computation of the free en¬ 
ergy minima, transition states and minimum energy pathways 
are described in the appendix. 


III. MINIMUM FREE ENERGY STATES 

We first study the minimum free energy morphologies of 
the system as the surface anchoring W is varied. This is sim¬ 
ilar to the work reported in where the authors use contin¬ 
uation methods to compute stable equilibria of this system as 
W is varied. Throughout this paper, we fix the value of the 
dimensionless elastic constant /^ez = 4 x 10“^ (for typical 
values of the elastic constant and Landau-de Gennes bulk po¬ 
tential parameters, this describes a micron-scale well) and use 
a square lattice with 150 x 150 grid points. 

In agreement with previous studies Halil , we find that 
for large anchoring strengths i.e. for W > 1.4 x 10“^, there 
are two classes of minima, namely diagonal and rotated states. 
We check the stability of the solutions by numerical computa¬ 
tions of the Hessian of the Landau-de Gennes energy and its 



FIG. 1: Local free energy minima for a square device for W > 
1.4 X 10“^. There are two classes of configurations: (a-b) diag¬ 
onal and (c-f) rotated states. Each double arrow indicates the two 
minima are connected by a single transition state, thus visualising 
the connectivity of the free energy landscape. 



FIG. 2: Minimum free energy diagonal configurations in a square 
device fovW < 1.4 x 10“^. The double arrow indicates that the 
two minima are connected by a single transition state. 


eigenvalues at a given solution, as is standard for numerical 
stability analysis. Due to the system symmetry, there are two 
equivalent diagonal states and four equivalent rotated states. 
These states are shown in Fig.[T] 

For weak anchoring, W < 1.4 x 10“^, the rotated states 
are no longer stable minima. This is consistent with the bi¬ 
furcation diagram in (1] where the authors do not find rotated 
solutions for small W. However, the rotated solutions do sur¬ 
vive as transition states in the free energy landscape for weak 
anchoring, connecting the stable diagonal solutions. This is 
further explained in the next section and the two stable min¬ 
ima (diagonal states) are shown in Fig.[2l 

Each double arrow in Figs. [T] and [2| indicates that there ex¬ 
ists at least one transition pathway between the connected 
minima and we discuss the transition pathways in the next 
section. 
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IV. TRANSITION PATHWAYS 


Previous studies have been limited to free energy minima 
in the planar bistable device HH. In this paper we system¬ 
atically compute the transition states, transition pathways and 
optimal transition pathways as a function of the surface an¬ 
choring parameter W. To the best of our knowledge, this has 
not been reported elsewhere in the literature. 

We distinguish three separate regimes depending on the val¬ 
ues of W : (i) strong anchoring regime, (ii) moderate anchor¬ 
ing regime and (iii) weak anchoring regime. 

Regime I: Strong anchoring regime - The first regime is 
the strong anchoring limit, which we find to occur for W > 
6.5 X 10“^. The optimal transition pathways are independent 
of W in this regime whilst they are sensitive to W for weaker 
anchoring. The computations in this regime are in fact very 
similar to the Dirichlet case considered in m where Q is 
fixed to be q on the square edges. 

We first make a few comments about the vertex defects in 
the diagonal and rotated solutions, both of which are local 
free energy minima in this regime. There are two types of 
point defects in these solutions: (i) splay defects wherein the 
director splays outward from the vertex and (ii) bend defects, 
for which the director bends between a pair of intersecting 
edges. 

For a diagonal solution, there are two diagonally opposite 
splay defects and two diagonally opposite bend defects and 
for a rotated solution, the two splay defects are connected by 
an edge, as are the two bend defects. For a rotated to diagonal 
transition, we need two defect transformations along an edge: 
a bend to splay defect transformation and a splay to bend de¬ 
fect transformation via director rotation. We obtain at least 
three different transition pathways between a rotated and a di¬ 
agonal solution, as illustrated in Figure [3] below, all of which 
are mediated by motion of point defects. 

The first mechanism, as shown in Fig.[3ta-c), is mediated by 
a — 1 /2 defect. The defect is created near the top left corner, 
propagates towards the bottom left corner and subsequently 
settles there to yield a bend defect. The numerical results sug¬ 
gest that the singular bend behaviour at the top left corner of 
the initial rotated solution is propelled into the square inte¬ 
rior leading to the creation of a —1/2 defect. As the —1/2 
defect propagates downwards, it rotates simultaneously and 
induces a global director rotation within the square domain. 
Finally, the point defect settles at the bottom left corner to 
become a bend defect and the system relaxes into a diagonal 
solution. The second mechanism. Fig. Od-f), is dominated 
by a H-1/2 defect created near the bottom left corner. The 
rotated solution has a splay defect at the bottom left corner; 
this defect is propelled into the square interior, it moves up¬ 
wards along the square edge and whilst moving, rotates and 
induces a global director rotation. It finally settles at the top 
left corner as a splay defect and the system relaxes into a di¬ 
agonal solution. These two mechanisms have identical free 
energy profiles, as shown in Figure 5. Further analysis using 
our numerical methods also shows that there can be other tran¬ 
sition pathways with higher free energy barriers. We interpret 
the optimal transition pathway as being the minimum energy 



FIG. 3: Possible transition pathways between rotated and diagonal 
states in the strong anchoring limit, W > 6.5x10“^. The transitions 
are mediated by a —1/2 defect, a +1/2 defect, and a pair of ±1/2 
defects respectively in panels (a-c), (d-f), and (g-i). 



FIG. 4: (color online) The order parameter plots, s{x,y), for the 
transition pathways shown in Fig. [3] s ^ 0 at the defect locations 
and s ^ 1 elsewhere. The transitions are mediated by a —1/2 de¬ 
fect, a +1 /2 defect, and a pair of ± 1 /2 defects respectively in panels 
(a-c), (d-f), and (g-i). 
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Normalised Path Length Normalised Path Length 


FIG. 5: (a) The free energy profiles for the three transition mech¬ 
anisms discussed in Fig. [3] The —1/2 and +1/2 defect mecha¬ 
nisms have degenerate free energy profiles, while the ±1/2 defect 
pair mechanism has a considerably higher free energy barrier, (b) 
The free energy profiles along the optimal transition pathway for dif¬ 
ferent values of the surface anchoring strength W. 


pathway with the smallest barrier and the least energetically 
expensive transition state. An example of a transition path¬ 
way with a higher free energy barrier is shown in Fig. [3](g-i), 
where the pathway is mediated by the creation and propaga¬ 
tion of a pair of ±1/2 defects. We also note that the transition 
pathways between other rotated and diagonal states can be re¬ 
lated by symmetry to the specific example shown in Fig.O In 
Fig. m we plot the order parameter along the transition path¬ 
ways, as shown in Fig. 0 The defects at the four corners, 
along with the bulk ±1/2 defects (mediating the transitions), 
are clearly identified by small order i.e. by small values for s 
such s = \/Qii + Qi 2 0 near the defects. 

The free energy profiles for the three transition mechanisms 
discussed in Fig. [3] are shown in Fig. [3a). The x-axis corre¬ 
sponds to the ''normalisedpath length'' between configuration 
m along the pathway and the initial rotated state, 


N 


1 1/2 




'ij \n 
a^) 


n=0 a (3 ij 


(5) 


The normalization constant N is chosen such that the total 
path length between the rotated and diagonal states is 1. For 
the first two mechanisms described above, the profiles are vir¬ 
tually indistinguishable and the maxima in the free energy 
profiles correspond to the saddle points (transition states) in 
the free energy landscape of the system. The transition state 
configurations are depicted in Figs. Ob), (e), and (h). It is 
worth noting that we do not find any direct pathway between 
two diagonal or two rotated states. A transition between two 
diagonal states or two rotated states is always composed of a 
sequence of rotated-diagonal transitions. 

Regime II: Moderate Anchoring - We label the range 1.4 x 
10-^ < W <6.5x10 ^as the medium anchoring regime. 
There are three key features of this regime: (i) the diagonal 
and rotated solutions survive as local free energy minima, (ii) 
the free energy barrier (along minimum energy pathways) de¬ 
creases monotonically with W and (iii) the optimal transition 
pathways do not feature defects. 

This regime is further divided into two sub-regimes (a) 
2.4 X 10-^ < VF < 6.5 X 10“^ and (b) 1.4 x 10“^ < VF < 



FIG. 6: Two competing transition pathways in the medium anchor¬ 
ing regime, 2.4 X 10“^ < W < 6.5 X 10~^. The molecules rotate 
clockwise in the first mechanism (a-c). The director rotation sweeps 
upward. In the second mechanism (d-f), the rotation is also clock¬ 
wise, but it starts from the top left corner and propagates downward. 


2.4 X 10“^. For case (a), the transition from a rotated to a 
diagonal solution is due to a "sequential" director rotation on 
one side of the square device, without any transient defects. 
Heuristically, it is energetically preferable to break the tan¬ 
gential anchoring along a square edge leading to director rota¬ 
tion along an edge as opposed to the creation of ±1/2 defects. 
As in the strong anchoring regime, we find that there are two 
mechanisms with degenerate free energy profiles (including 
the free energy barrier). 

In the first mechanism, the director starts to rotate clock¬ 
wise from the bottom left corner, breaking the tangential an¬ 
choring along the left vertical edge. The rotation then propa¬ 
gates upwards, as shown in Fig. Oa-c) and this rotation suf¬ 
fices to transform the splay defect at the bottom left corner 
into a bend defect and conversely, the bend defect at the top 
left corner into a splay defect. The final state is a diagonal so¬ 
lution. The second mechanism is shown in Fig. (Sd-f). The 
director rotates clockwise from the top left corner and the 
sweeping motion propagates downwards, inducing a global 
director rotation and a transition from a rotated solution to a 
diagonal solution. The configurations in panels (b) and (e) 
correspond to the transition states and the two mechanisms 
have identical free energy profiles. 

For weaker anchoring, for 1.4 x 10“^ < FF < 2.4 x 10“^, 
the penalty for breaking tangential anchoring is weak and the 
rotated to diagonal transition follows from a "global" director 
rotation along the left vertical edge. As shown in Fig.Oa-d), 
the director rotates simultaneously along the entire length of 
the left vertical edge and this rotation induces a global director 
distortion within the square, leading to a rotated to diagonal 
transition. 

Regime III: Weak Anchoring - As the surface anchoring 
strength W further decreases, for W < 1.4 x 10“^, the ro¬ 
tated solutions are no longer stable. The disappearance of 
the rotated states correspond to cusp catastrophes ill , where 
typically a minimum and two transition states coalesce to a 
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FIG. 7: The minimum energy transition pathway between a rotated and a diagonal state in the medium anchoring regime, 1.4x10 ^ <W < 
2.4 X 10“^. Panel (b) corresponds to the transition state along this pathway. 


single point. This point, reminiscent of the rotated state (as 
such, there are four of them), is a saddle point (transition 
state) between the two stable diagonal solutions. The transi¬ 
tion state in panel Fig.[8](c) resembles a rotated solution, with 
a uniform director profile in the square interior accompanied 
by two transition layers near the vertical edges. The corre¬ 
sponding free energy profile can be seen from Fig. [3b) for 
W = 1.0 X 10“^. In this regime, the minimum energy path¬ 
ways may therefore be represented by the diagram in Fig. |2| 
where the double arrow represents four equivalent transition 
pathways. 

The movies for all the possible transition pathways de¬ 
scribed above are available as supporting information. 


V. CONCLUSIONS 


We have systematically studied the free energy landscape 
of a multistable nematic liquid crystal device as we vary the 
strength of the surface anchoring potential. Previous work fo¬ 
cuses on the free energy minima and we take this work further 
by computing transition states, transition pathways and opti¬ 
mal transition pathways between minima, along with the free 
energy barriers. 

We have classified the system behaviour into three cate¬ 
gories, according to surface anchoring strength and qualita¬ 
tive properties of transition pathways. In the strong anchoring 
regime, W > 6.5 x 10“^, we obtain two distinct classes of 
free energy minima: diagonal and rotated solutions, as ex¬ 
pected from previous work. The optimal transition pathways 
and the corresponding transition states feature either a +1/2 
or a —1/2 defect, localized near an edge. It is worth noting 
that these transition states bear strong resemblance to exper¬ 
imental observations of states with internal defects reported 
in |@], suggesting possible metastability of transition states. 
Further, there are multiple transition pathways between free 
energy minima, which can be mediated by multiple defects but 
with higher free energy barriers than the optimal pathways. 

The existence of a transition state or a saddle point-type 
critical point of the Landau-de Gennes energy (jU), connecting 
the stable diagonal and rotated solutions in the strong anchor¬ 
ing regime, can be theoretically justified from the celebrated 
Mountain Pass Theorem H. In the limit of infinite anchor¬ 
ing, one can define a topological classification scheme for 
two-dimensional director fields (see S for a similar three- 
dimensional topological classification), for which the diago¬ 


nal and rotated solutions belong to different topological sec¬ 
tors. This can also be seen more intuitively from the arrange¬ 
ment of the splay and bend point defects in the diagonal and 
rotated solutions respectively. The Mountain Pass Theorem 
guarantees the existence of at least one critical point of the 
Landau-de Gennes energy connecting the topologically dis¬ 
tinct diagonal and rotated solutions. In our framework. Moun¬ 
tain Pass solutions correspond to transition states along opti¬ 
mal transition pathways. 

In the moderate anchoring regime, for W <6.5x10“^, we 
do not observe any defect along the optimal transition path¬ 
ways. The rotated-to-diagonal transition is driven by clock¬ 
wise director rotation (and anti-clockwise for the opposite 
diagonal-to-rotated transition) localized along an edge. For 
2.4 X 10“^ < FF < 6.5 X 10“^, the rotation starts from one 
corner and propagates to another corner, along an edge. For 
weaker anchoring strengths (W < 2.4 x 0“^), it is preferable 
for the director to rotate along an entire edge, since the energy 
penalty for violating tangent anchoring is relatively weak. 

For W < 1.4 X 10 we capture the cusp catastrophe 
events as the rotated states become unstable and act as tran¬ 
sition states for the transition pathways connecting two stable 
diagonal solutions. This was previously not noted in the liter¬ 
ature. 


We note that the same effect i.e. loss of stability of the ro¬ 
tated solutions, can also be achieved by continuously decreas¬ 
ing the square size whilst maintaining strong anchoring on the 
square edges. In ifT^ . the authors study the critical points 
of the two-dimensional Landau-de Gennes energy as a func¬ 
tion of the square size, with strong anchoring on the edges. 
For micron-sized wells, they recover the diagonal and rotated 
solutions as we do. For a critical well size, the rotated so¬ 
lutions lose stability (but exist as non-minimizing solutions) 
whilst the diagonal solutions retain stability. For nano-scale 
wells, there is a unique critical point, referred to as the order- 
reconstruction solution, by analogy with similar findings in 


For moderate and weak anchoring, the height of the free 
energy barrier decreases monotonically with W whilst the 
height of the free energy barrier is independent of W in the 
strong anchoring regime. Secondly, the configurations of the 
transition states and the free energy barriers depend strongly 
on the anchoring strength, for moderate and weak anchoring. 
These results can therefore be exploited to fine-tune the de¬ 
sirable surface properties i.e. use an anchoring strength that 
yields the optimal free energy barrier for device applications. 
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FIG. 8: One of four equivalent optimal transition pathways between two diagonal states in the very weak anchoring regime, VF < 1.4 x 10 ^ 
Panel (c) corresponds to the transition state, which is reminiscent of a rotated state. 


For example, a very low barrier compromises the stability of 
the minimum free energy configurations and the device per¬ 
formance. A very high barrier is an impedance to realistic 
switching, necessitating large power input or strong external 
electric fields. 

We will now discuss how the dimensionless free energy can 
be converted to SI unit (Joules). Note that the free energy de¬ 
fined in Eq.(|4|) is in two dimensions (i.e. free energy per unit 
thickness), while all practical experiments are, of course, in 
three dimensions. As such, we need to account for the thick¬ 
ness of the device, A 2 ;, such that ^ = L‘^/S.Z/C)^, where 
^ is the dimensionless free energy defined in Eq. 0. The 
parameters a and C are experimentally measurable quanti¬ 
ties. As an example, for the common liquid crystal mate¬ 
rial MBBA the characteristic bulk constants are of order 
a ^ 4.2 X 10^ J/m^ and C ^ 3.5 x 10^ J/m^. To estimate a, 
we have assumed that the temperature is 1° below the critical 
nematic transition temperature 46°C ). Erom Eig. [5l the 
typical free energy barrier in dimensionless units is of order 
10“^. Assuming the thickness of the device is 1 micron and 
the area of the device is (100 micron)^, the typical free energy 
barrier in SI unit is therefore of order 5.1 x 10“^^ J. This is, of 
course, much larger than ^ 4.4 x 10“^^ J. On the other 
hand, if the device dimension can be miniaturized to 1 micron 
X 1 micron x 100 nm, then the typical free energy barrier be¬ 
comes of order 5.1 x 10“^^ J, comparable to /c^T. The free 
energy barriers for other materials and device dimensions can 
be computed in an analogous way. 

Our methods can be extended to three dimensions and to 
include external fields. Eor example, we could use similar 
methods to study transition pathways in three-dimensional 
set-ups, as studied in |@], where the authors employ a three- 
dimensional Landau-de Gennes modelling approach and find 
a new biaxial order-reconstruction pattern for shallow nano¬ 
scale rectangular wells with strong anchoring. We speculate 
that the biaxial order-reconstruction pattern loses stability for 
larger wells but exists as a transition state connecting diagonal 
and rotated solutions for large micron-scale wells with strong 
anchoring. 

Similarly, we can use these numerical methods to describe 
more complex geometries and boundary conditions, for exam¬ 
ple the Zenithally Bistable Device studied in ifTsIl . They are 
also compatible with other simulation techniques tailored to 
the dynamics of liquid crystalline systems or generally com¬ 
plex fiuids problems (such as finite element, lattice Boltzmann 
method, etc). These numerical methods promise to be a valu¬ 


able design or optimization tool, allowing us to obtain a de¬ 
tailed picture of the free energy pathways of a complex sys¬ 
tem, and subsequently how the pathways evolves as the sys¬ 
tem parameters are varied. 


VI. APPENDIX: COMPUTATIONAL TECHNIQUES 

Here we summarise the most salient features of the compu¬ 
tational techniques. A more detailed description and imple¬ 
mentation of the methods can be found in ifo . 


A. Discretization of the Free Energy Functional 

We discretize the Landau-de Gennes free energy on a 
square grid. Specific stencils are needed to approximate the 
derivatives and it is important that they are at least second or¬ 
der accurate. We have used 


|VQ*^' 1 2 




+ -Q%f 

2(Ax)^ 

+ -Q%? 

2(A2/)2 


where the superscripts i and j label the lattice points in two 
dimensions. We have used a 150x150 grid based on the results 
of Luo et al. (H], where they carried out convergence tests 
on the energy-minimizing solutions in the Landau-de Gennes 
framework. 

Minimum Free Energy States - The first step in a survey 
of the free energy landscape is to compute the majority, if 
not all, of the possible minima in the system. We follow a 
stochastic approach, where each step consists of a trial move 
followed by an energy minimization. It is identical to a basin¬ 
hopping algorithm 111, [H] at infinite temperature, such that 
every minimum state found is recorded. The simplest trial 
move consists of random perturbations of the lattice field val¬ 
ues, (QOa/3 ~ ^ is a random number be¬ 

tween -1.0 and 1.0, and 0 is the amplitude of the perturba¬ 
tion, which we usually take to be 0 = 0.5. The energy min¬ 
imization is carried out using the limited-memory Broyden- 
Eletcher-Goldfarb-Shanno (LBEGS) algorithm llll [U. In 
this paper, we typically take 500 (basin-hopping) steps to sam¬ 
ple the minima of the system. 
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B. Transition States and Minimum Energy Pathways 

Given the minimum free energy states, we numerically 
compute the transition states (saddle points in the free energy 
landscapes) using a combination of the doubly-nudged elastic 
band (DNEB) meth od [1^ and the hybrid eigenvector¬ 
following technique 11^ . Transition states are special saddle 
points or critical points in the free energy landscapes where 
the energy gradients are zero in all eigendirections and one 
of the eigenvalues of the Hessian is negative. This is to be 
contrasted to minima, where the gradients are zero and the 
eigenvalues are positive in all eigendirections. 

The DNEB method is a double-ended search algorithm for 
finding transition states and minimum energy pathways be¬ 
tween any two pair of minima. Minimum energy pathways 
have an important feature such that every point on such a path¬ 
way is a minimum in all but a single distinguished direction 
in the phase space. As a consequence, a maximum along the 
minimum energy pathway corresponds to a transition state. 

InDNEB, a set of images, {E^, ..., F^}, areplaced be¬ 
tween the two endpoints (minima), F^ and The sym¬ 

bol r represents all the degrees of freedom in the system, i.e. 
{Q^/ 3 }. We typically use 30 images, and they are initialised 
by taking a linear interpolation between the two endpoints. 

The images are relaxed using Landau-de Gennes energy 
gradients that include contributions from two components. 
The first contribution, g, comes from the derivatives of the 
free energy functional with respect to the lattice degrees of 
freedom, The second contribution is a spring 

force, g, which is required to keep the images roughly equidis¬ 
tant. The following spring potential is applied on every image 

^spring = f ((«“’-)' - (s“’ + )2) . (6) 

and are respectively the distances to the left and 
right images, ( 5 ^’^)^ = — F^p. The coefficient k is 

the spring constant. 

Using the true and spring gradients, without any projection. 


often results in two issues (i) corner-cutting, where 

images are pulled away from the minimum energy path; and 
(ii) sliding-down problems, where images slide down from 
barrier regions. In DNEB, we use two projections. Eirst, we 
only retain the components of the true gradient that are per¬ 
pendicular to the unit tangent vector f ^, 

gj =g"-(g"-r")f-. (7) 

Secondly, the following component of the spring constant gra¬ 
dient is retained, 

gDNEB = + g± - (Sl ■ g“)gi, (8) 

where g? = g*^ — (g*^ • t^)t^. The unit tangent f is defined 
as in 11^ . As a convergence criterion, we stop the DNEB 
run when the root mean square of the gradients as defined 
in Eqs. © and ^ in the manuscript is < 10 which in 
our experience is more than enough to ensure the calculations 
have converged. 

A maximum in the DNEB path corresponds to a transition 
state candidate, which we then refine by applying the hybrid 
eigenvector-following technique. We use a Rayleigh-Ritz ap¬ 
proach im, based on gradients of the Landau-de Gennes en¬ 
ergy, to compute the transition states since this method is more 
efficient for large systems. 

Once the transition states are found, small displacements 
are applied in the two downhill directions. Energy min¬ 
imizations are again carried out using the LBEGS algo¬ 
rithm imili. Given a minima - transition state - minima 
triplet, this run yields the transition pathway for this triplet 
and the corresponding free energy barrier. We check that this 
pathway is consistent with the one obtained from DNEB cal¬ 
culations. The free energy barriers could be interpreted as the 
free energy differences between the transition state and the 
free energy minima. There are typically multiple transition 
pathways between a pair of free energy minima. The optimal 
pathway is the pathway with the smallest free energy barrier. 
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